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Exact equations are derived that relate velocity structure functions of arbitrary order with other 
statistics. "Exact" means that no approximation is used except that the Navier-Stokes equation and 
incompressibility condition are assumed to be accurate. The exact equations are used to determine 
the structure-function equations of all orders for locally homogeneous but anisotropic turbulence as 
well as for the locally isotropic case. These equations can be used for investigating the approach to 
local homogeneity and to local isotropy as well as the balance of the equations and identification of 
scaling ranges. 



I. INTRODUCTION 



Full mathematical exposition on the topic of structure-function equations is given here. A brief summary 
of results derived here will appear in the Journal of Fluid Mechanics in the paper "Equations relating structure 
functions of all orders." The two sections below are sufficiently similar to that paper so as to guide the reader to 
the relevant mathematical details, much of which resides in the appendices herein. The two sections below contain 
more mathematical detail than does that paper. The derivation of the structure function equations of all orders 
produces substantial mathematical detail. This is true for reduction of the viscous term and for the term involving 
the pressure gradient when deriving the exact equations. Applying isotropic formulas for structure functions of 
arbitrary order requires the invention of new notation and much use of combinatorial analysis. The divergence 
and Laplacian operating on isotropic formulas necessarily appear in the equations; evaluation of which requires the 
derivation of many identities. Finally, matrix-based algorithms are invented such that the isotropic formulas for the 
divergence and Laplacian of isotropic tensors of any order can be generated by computer. 

There is some difference in notation between the paper and this document. In the paper, a component 
of a structure function is denoted by D^i^^ jsi^ i^^-^, whereas here it is denoted by the more complicated notation 
D[N:Ni,N2,N3]- The reason for the more complicated notation here is to avoid ambiguity at several places in the mathe- 
matics. In the paper, the components of the tensor { W[7v_2P] (r) S[2P] } are denoted by { W[jv-2P] (r) S[2P] } jy n ]' 

Here, there is no symbolic distinction between the tensor {W[n-2P] (r) (^pp]} and its components. The distinction 
is implied by the context. 

II. EXACT TWO-POINT EQUATIONS 

The Navier-Stokes equation for velocity component Ui(x, t) and the incompressibility condition are 

dtUiiyi, t) + u„(x, t)dx^Ui{yi, t) = -dx,p{'x,t) + vdx„dx„Ui{y.,t) , and 5a;„u„(x, t) = 0, (1) 

where p(x, t) is the pressure divided by the density (density is constant), v is kinematic viscosity, and d denotes partial 
differentiation with respect to its subscript variable. Summation is implied by repeated Roman indexes. Consider 
another point x' such that x' and x are independent variables. For brevity, let Ui = Ui{ii.^t), u[ = Ui{x',t), etc. 
Require that x and x' have no relative motion. Then dx^Uj = 0, dx'Uj — 0, etc., and dt is performed with both held 
X and x' fixed. Subtracting (|^) at x' from (|^) at x and using the aforementioned properties gives 

dtVi + Undx^Vi + u'j^dx'^Vi = -Pi + V (dx„dx„Vi + dx'^dx'^Vi) , (2) 
where Vi = Ui - u'^ , P^ = {dx,p - dx'.p'^ ■ (3) 

Change independent variables from x and x' to the sum and difference independent variables: 

X = (x + x') /2 and r = x — x', and define r = |r| . (4) 
The relationship between the partial derivatives is 
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dx, = dr, + ^dxi , dxi = -dr, + ^dx, , dx, = dxi + dx'^ , = ^ [dx, - dx'^ ■ (5) 

The change of variables organizes the equations in a revealing way because of the following properties. In the case of 
homogeneous turbulence, dxi operating on a statistic produces zero because that derivative is the rate of change with 
respect to the place where the measurement is performed. Consider a term in an equation composed of dxi operating 
on a statistic. For locally homogeneous turbulence, that term becomes negligible as r is decreased relative to the 
integral scale. For the homogeneous and locally homogeneous cases, the statistical equations retain their dependence 
on r, which is the displacement vector of two points of measurement. Using (||), (||) becomes 

dtVi + Undx„Vi + Vndr,Vi -Pi + V {dx„dx„Vi + dx'^dx'^Vi) , (6) 

where ?7„ = [ui + u[) /2. (7) 

Now multiply (|^) by the product VjVk ■ ■ ■ vi, which contains — 1 factors of velocity difference, each factor 
having a distinct index. Sum the N such equations as required to produce symmetry under interchange of each pair 
of indexes, excluding the summation index n. French braces, i.e., {o}, denote the sum of all terms of a given type 
that produce symmetry under interchange of each pair of indexes. The differentiation chain rule gives 

{vjVk ■ ■ ■ vidtVi} = dt (vjVk ■ ■ ■ viVi) , (8) 

{VjVk ■ ■ ■ VlUndx„Vi} = Undx„ (VjVk ■ ■ ■ VlVi) = dx„ {UnVjVk ■ ■ ■ VlVi) , (9) 
{VjVk ■ ■ ■ VlV„dr„Vi} = Vndr„ {VjVk ' ' ' VlVi) = 5r„ {VnVjVk ■ ■ ■ VlVi) . (10) 

The right-most expressions in (^) and (|lO|) follow from the incompressibility property obtained from (|^) and the fact 
that dxiU'j = 0, dx'Uj = 0, namely 

dx„Un = 0, dx„Vn = 0, dr^Un = 0, 9r„Wn = 0. (11) 

The viscous term in (|^) produces i> {vjVk ■ • ■ vi [dx„dx„Vi -\- dx'^dx'^Vi)Y, this expression is treated in Appendix A. 
These results give 



dt [Vj • • • Wj) + dx„ {UnVj ■ ■ - Vi) + dr„ {VnVj ■ ■ ■ Vi) 



{v, ■ ■ ■ viP,} + 2v 



dr„dr„ + \dx„dx„ ) [vj ■ ■ ■ Vi) - {vk ■ ■ ■ vie^j) 



(12) 
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where = {dx^u^) (dx^uj) + {dx'^u'i) {dx'^u'j) = {dxr,'"i) (dx^Vj) + [dx'^Vi) {dx'^Vj) . (13) 

The quantity {vj ■ ■ ■ viPi} can be expressed differently on the basis that (^ allows Pi to be written as Pi = 
dxi [p — p')- The derivation is in Appendix B; the alternative formula is 

{vjVk ■ ■ ■ viP,} = {dx^ [vjVk ---viip- p')]} -{N-l){p- p') {{s,j - s'y) Vfe • • • V/} , (14) 

where the rate of strain tensor Sij is defined by 

Sij = {dx.Uj + dxjUi) /2. (15) 



III. AVERAGED EQUATIONS 

Consider the ensemble average because it commutes with temporal and spatial derivatives. The above notation 
of explicit indexes is burdensome. Because the tensors are symmetric, it suffices to show only the number of indexes. 
Define the following statistical tensors, which are symmetric under interchange of any pair of indexes, excluding the 
summation index n in the definition of Fjyv^]^]: 

D[7V] = {Vj ■ ■ - V^) , F[N+1] = {UnVj ■ ■ ■ Vi) , T[N] = {{Vj ■ ■ ■ ViPi}) , E[Ar] EE {{vk ■ ■ ■ Vidj}) , (16) 

where angle brackets () denote the ensemble average, and the subscripts N and + 1 within square brackets denote 
the number of indexes. The argument list (X, r, t) is understood for each tensor. The left-hand sides of each 
definition in (|l^) are in implicit-index notation for which only the number of indexes is given; the right-hand sides in 
( p^ ) are in explicit-index notation. The ensemble average of (|l^) is 



dt'D^i^] + Vx • F[Ar+i] + Vr • D[^_|_i] — -T[Ar] + 2j/ 



4 ■ X ] ^[N] - -C^IW] 



(17) 



where, Vx • F[Ar+i] = 9x„ ^nVj ■ ■ -v,) ,Vr ■ D[Ar+i] = dr„ {vnVj ■ ■ ■ Vi) ,Vl = i9r„i9r„, Vx = dx„dx„- The notations 
Vx-, Vx, Vr', and V^ are the divergence and Laplacian operators in X-space and r-space, respectively. 
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A. HOMOGENEOUS AND LOCALLY HOMOGENEOUS TURBULENCE 



Consider homogeneous turbulence and locally homogeneous turbulence; the latter applies for small r and large 
Reynolds number. The variation of the statistics with the location of measurement or of evaluation is neglected 
for these cases. That location being X, the result of Vx' operating on a statistic is neglected. Thus the terms 
Vx • F[jv_|_i] and jVxDfAr] are neglected in (p7|); then ( |l7| ) becomes 

dtB^N] + Vr • D[w+i] = -T[Ar] + 2iy [wIB^n] - ^n]] ■ (18) 

Because the X-dependence is neglected, the argument list (r, t) is understood for each tensor. The ensemble av- 
erage of (|lj) contains {dxi [{vjVk ■ ■ ■ vi {p — p')}]), which can be written as the sum of TV — 1 statistics of the form 
{jvjVk ■ ■ ■ vi {p ~ p')}) operated upon by the X-space gradient. Since such X-space derivative terms are neglected, 
(|lj) gives the alternative that 

T[^] = - (TV - 1) ((p -p') {(%■ - 4) • ■ • ^'}> ■ (19) 

Locally homogeneous turbulence is also locally stationary such that the term c^tDj^v] in ( p^ may be neglected. 
However, dt'D[j\j] is not necessarily negligible for homogeneous turbulence. 



B. ISOTROPIC AND LOCALLY ISOTROPIC TURBULENCE 



Consider isotropic turbulence and locally isotropic turbulence; the latter applies for small r and large Reynolds 
number. The tensors D[7v], T[jv], and E[Ar] in (^6|) obey the isotropic formula. The Kronecker delta Sij = 1 if 
i = j and Sij — ii i ^ j. Let S[2P] denote the product of P Kronecker deltas having 2P distinct indexes, and let 
W[7v] (r) denote the product of N factors ^ each with a distinct index; the argument r is omitted when clarity does 
not suffer. Because each tensor in ( |l6|) is symmetric under interchange of any two indexes, their isotropic formulas 
are particularly simple. Each formula is a the sum of M + 1 terms where 

M = N/2 if N is even, and Af (iV - 1) /2 if N is odd. (20) 

Each term is the product of a distinct scalar function with a W[jv] and a S^]- From one term to the next, a pair 
of indexes is transferred from a W[jv] to a (5[2P]; examples are given in ( j65| - |67[ ) of Appendix E. For the tensor D^jv], 
denote the Pth scalar function by Dn,p (r, t). Thus the scalar functions belonging to the isotropic formulas for T[jv] , 
E[^], and D[jv+i] are denoted by T^y.p {r, t), Ejq p (r, t), and _Djv+i.p {^i t), respectively. The scalar functions depend 
on the magnitude of the spacing r rather than on the vector spacing r. The isotropic formula for D[jv] is 

M 

D[jv] (r, t)=Yl P { W[Ar-2P] (r) ^pp] } , (21) 

P=0 

and the isotropic formulas for T[n] and E[^] have the analogous notation. Recall that {o} denotes the sum of all 
terms of a given type that produce symmetry under interchange of each pair of indexes. Henceforth, the argument 
list (r, t) will be deleted. 

A special Cartesian coordinate system is typically used because it simplifies the isotropic formula. This 
coordinate system has the positive 1-axis parallel to the direction of r, and the 2- and 3-axes are therefore perpendicular 
to r. Let A^i, N2, and be the number of indexes of a component of D[Ar] that are 1, 2, and 3, respectively; such 
that N = Ni + N2 + N3. Because of symmetry, the order of indexes is immaterial such that a component of Dj^r] 
can be identified by A'^i, and N3. Thus, denote a component of D[Ar] by D^p^.j^-^ pf^j^^^, which is a function 
of r and t. The projection of ( pT] ) using Ni, N2, and N3 unit vectors in the directions of the 1-, 2-, and 3-axes, 
respectively, results in the component D^j^j.j^-^ j^^ pf^] on the left-hand side of (pT|), and numerical values of the projection 
of {W[^_2P] (r)<5[2P]} appear on the right-hand side. Henceforth the word "projection" will be omitted for brevity. 
Those values of the coefficients {W[jv-2P] (r) ^[2P]} in (21) are needed; the values obtained for the special coordinate 



system are determined in Appendix C; they are, from (44 



if 2P < iV2 -I- N3 then { W[jv-2P]^[2P] } = 0; otherwise, (22) 



{W[Ar_2P]'5[2P]} = A^l!iV2!^3!/ 



(^-2P)!2-(f ! f ! P-f )! 



(23) 
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By applying ( pT| ) and (]2^p3|) for all combinations of indexes, one can determine which components -DjTViTVi.iVa.Wa] 
are zero and which are nonzero, identify M + 1 linearly independent equations that determine the D^ p in terms 
of M + 1 of the -D[7V:7Vi,Af2.Af3]' S^id algebraic relationships between the remaining nonzero D[N,N^^ff^^ff^]. The 
derivations are in Appendix D; a summary follows. 

A component Di^p^.j^-^ pf^ pf^] is nonzero only if both N2 and N3 are even and when iVi is odd if N is odd and 
when A^i is even if N is even. Thereby, (M + 1) {M + 2) /2 components are nonzero. There are 3^ components of 
D[jv]; thus the other 3^ — (M + 1) (M + 2) /2 components are zero. 

There exists exactly (M + 1) M/2 kinematic relationships among the nonzero components of D[jv]. For each 
of the M + 1 cases of iVi , these relationships are expressed by the proportionality 
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[N:Ni,2L,Q] 



D 



[N:Ni,2L~2,2] ■ £'[jV:Ari,2L~4,4] 
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[N:Ni,0,2L] 



[(2L)!0!/L!0!] : [{2L - 2)!2!/ {L ~ 1)!1!] : [(2L - 4)!4!/ {L - 2)!2!] 



: [0! (2L)!/0!i 



(24) 



Previously, only one such kinematic relationship was known (Millionshtchikov 1941). For = 4, ( |24|) gives I?[4:o.4,o] ■ 
^[4:0.2,2] '■ D[4.Q Q 4] = 12 : 4 : 12. In explicit-index notation this can be written as 132222 = 31)2233 = £'3333, which 
was discovered by Millionshtchikov (1941). Now, all such relationships are known. 

There remain M + 1 linearly independent nonzero components of D[^]. This must be so because there are 
M + 1 terms in (|2l|), and the M + 1 scalar functions Dj^ p therein must be related to M + 1 components. Consider 



the M + 1 linearly independent equations that determine the Dn^p in terms of M + 1 of the D 



simplicity, the chosen components can all have N^^ = 0, i.e. 



[N:Ni,N2,N3\- 



For 



2,2,0] I ^[Af:Af-4,4,0]; 



D 



[N:N-2M,2M,0\- 



the choice of linearly independent components can be 
As described above, projections of ( pT[ ) result in the 



chosen components on the left-hand side and algebraic equations on the right-hand side. These equations can be 
expressed in matrix form and solved by matrix inversion methods; the result is given in ( p8[ ) of Appendix F. Given 
experimental or DNS data or a theoretical formula for the chosen components, the solution in (g8|) determines the 
functions -Dat.p in (|2^); then ( pl| ) completely specifies the tensor Dj^v]- The matrix algorithm is an efficient means 
of determining isotropic expressions for the terms Vr • D^^v-i-i] and V^D[7v] in (p^. Those algorithms are given in 
Appendix F. From the example for A = 2 in Appendix F, use of the matrix algorithm and the isotropic formulas in 
(Es) gives the two scalar equations 



dtDii +[dr + -] Dm - -D122 
r r 



-Til + 2v 



2v 



dtD22 + I 9r + - I £'122 = -T22 + 2v 
r 



2v 



2 4 \ 4 

dl + -dr - — ]Dii + -^D22 - Ell 



dlDii + -drDii + — {D22 - Dii] 



4£/3, 



^Dii + (92 + _ ^ I - £;22 



dlD22 + -drD22 - \ {D22 - Dii) 



4e/3, 



(25) 



(26) 



where use was made of the fact (Hill, 1997) that local isotropy gives Tu ~ T22 = and 2uEii — 2vE22 = 4£/3 where 
e is the average energy dissipation rate per unit mass of fluid. Now, ( psfp^ ) are the same as equations (43-44) of 
Hill (1997), and Hill (1997) shows how these equations lead to Kolmogorov's equation and his 4/5 law. From the 
example for A^ = 3 in Appendix F, 



dtDiii 
dtD 122 4 



dr 



dr + - ) i^iiii - -i:'ii22 = -Till + 2iy[C- Em] . 
r J r 

4\ 4 

- £'1122 — 7;— £'2222 — —2^122 + 2l'[B — i?122 

r / Sr 



(27) 
(28) 



where 



-dr + df) Dm, and B = i f 4 - + 5d^ + rdl I Dm 



The incompressibility condition, D122 = J_(-Diii +r9r-Diii), was substituted in ( |113| ) to obtain (|2 
algorithm is checked by the fact that (^-29) are the same as given by Hill and Boratav (2001). 
The equations for A^ = 4 are 



(29) 

The matrix 



4 



dtDnn + { 9r + ^ ] £'11111 



-£>ii 



122 



= -Tiin + 21/ 



14. 



10 



+ -dr 5- -Dun + -5--D1122 + 77-5--D2222 



-2z/Eiiii, 



(30) 



dtDii22 + ( + - 1 -D11122 — ^-Dl2222 = 



1122 



2i/ 



2 

^i?iiii 



52 ^^2^2 



34 ■ 

-D1122 + ^^-D2222 



1122, 



(31) 



/ 6\ 

9tD2222 + [^^~^~) -^12222 = —22222 + 2l/ 



3r 



2222 



- 2i'E: 



2222- 



(32) 



Since these equations have a repetitive structure, it suffices to give the divergence and Laplacian terms. For iV = 5 
to 7 these are, respectively: 




{dr + 7) Dm 111 — ^-Diii 122 
{dr + 7) -Dill 122 ~ 7-D112 222 

{dr + f ) -D112222 — ■^-D222 222 



2u 



{d'r 



1 111 

_6 
5r 



la, 
+ 

-Dii 122 



-Dlllll — 7I-D11122 + |t-Di2 222 



10 



^9. i?i 



94 

5r^ 



12222 



+ (-^ + 52 + ^5,) i?i2 222 



111 122 



-Dun 111 — 

-Dun 122 ~ ||-Diii2 222 
-D1112222 — 57-D1222222 



2i/ 



V {dr + f ) -D1222 222 

/ (a2 + fa.-li)Dnini 
^Z)nini + (-l^ + 9.^ ' 



5r- 



■D 



br^ 
111 122 



r 

Idr) 



111 122 - 
-Dill 122 



+ (- 



/ {dr + 7) -Dn 111 111 - 
dr + 7) -Dn 111 122 - 
dr + I ) -Dn 112 222 - 

V {dr + f ) I?ll 222 222 - 

14^ 



112 
T5r^ 



dl 



920 

3; 

" T5r 



24 



12 222 



12 222 



14 ' 

^^Dlii 111 122 \ 



3f2 



416 

r-Da 

-"222 222 



ir;„2 -D222222 

248 , , 
12 222 — 75f5--i^222 222 

^ 4 £>o 



3;3--t-^112222 1572 -|- + r^r) -D222222 



11 112 222 



2!/ 



V 



r-Diin 111 + {—^^ + d. 
7^-1111122 + ( 15, 



22 222 222 

316 : 
r-^ 
^dr 



I {d1 + \dr — M) -Dim 111 — ^-Dnn 122 + ^^i^r-Diii2222 — 

21r 



111 122 



63r' 



D 



1112222 



' Dmi 122 +{-^+dl + Idr) I?1112 222 + 



1376 n 

-crj--t>'l222 222 

304 n 
" 15p--t>'l222 222 

1222 222 



7^2-01112 222 + (— sifr + 52 + ^dr) -D1222 222 
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V. APPENDIX A: THE VISCOUS TERM 



The quantity {vjVk • • ■ vidx„dx„Vi} requires special attention. Consider the repeated apphcation of the identity 
dx„dx^ ifg) = fdx„dx„g + gdx„dx,J + 2 {dx^f) {dx„g) (33) 

to the quantity 

dx„dx„ {VjVkVm ■ ■ - Vi) (34) 

for N factors of velocity difference in [vjVkVm ■ ■ ■ Vi). For the first application of ( ^3|) let / = Vj and let g be the 
remaining factors (wfcWm ■ • ■ Vi); this gives 

dx„dx„ (VjVkVm ■ ■ - Vi) ^ Vjdx„dx„ {VkVm ' -Vi) + [vkVm ' -Vi) dx„dx„Vj + 2 [dx„Vj] [dx„ {VkVm ' -Vi)] ■ (35) 

From the differentiation chain rule, dx^ {v^Vm ■ -Vi) is the sum of — 1 terms of the form (u,„ • -Vpdx^Vi). Thus, 
the right-most term in ( ^5|) is — 1 terms of the form 2vm ■ ■ ■ Vp {dx„Vi) (dx^Vj) each term containing N factors; 
two of those factors are distinguished by being derivatives of velocity differences. The second application of ( ^ ) is 
performed on Vjdx„dx„ (vkVm • 'Vi) in (p9), for which purpose f = Vk and g = {vm ■ -Vi); this gives 

Vjdx^dx„ {vkVm ■ ■ - Vi) ^ VjVkdx„dx„ (Wm ' -Vi) + Vj (w„ • -Vi) dx„dx„Vk + 2Vj [dx„Vk] [dx„ {Vm ' -Vt)] . 

The right-most term gives N — 2 terms of the form 2vjVm ■ ■ ■ Vp {dx„Vi) (dx^Vk) each term containing N factors. 

There are — 1 steps to complete reduction of the formula. The number of terms of the form 2vjVm ■ ■ ■ 
Vp (dx^Vi) {dx„Vk) is {N — 1) from the first step, (TV — 2) from the second step, etc. such that the total number of terms 
is (A^ - 1) -I- (A^ - 2) + • • • + (A^ - (A^ - 1)) = iV (iV - 1) /2. Now, A^ (TV - 1) /2 = (^) is the binomial coefficient equal 
to the number of ways of choosing two indexes from a set of A^ indexes; the quantities {dx^Vi) and (dx^Vj) in 2w,„ • 
• • Vp (dx„Vi) {dx„Vj) contain the chosen two indexes i and j. The (^) terms constitute 2 ■ ■ ■ Vp {dx„Vi) {dx„Vj)}. 
Because two factors of the form (vjVm ■ 'Vi) dx„dx„Vk appear in the last step, the total number of terms of the form 
(vj ■ ■ ■ Vn) dx„dx„Vi is A^. Not surprisingly, these N terms constitute {{vj ■ ■ ■ u„) dx„dx„Vi}, and A^ = (^) is the 
binomial coefficient equal to the number of ways of choosing one index from a set of A^ indexes, the quantity dx,^dx„Vi 
contains the chosen one index i. 
That is, for any A^ 

dx„dx„ (vj ■ ■ ■ v.,) = {{vj ■ ■ - vi) dx„dx„Vi} + 2{vk ■ ■ - vi {dx„Vi) (dx^Vj)} . (36) 

The left-hand side is symmetric under interchange of any pair of indexes (not including n because summation is 
implied over n), and the French brackets make the right-hand side likewise symmetric. 

Use of ( 13^ ) within the viscous term v {vjVk ■ ■ ■ vi [dx„dx„Vi + dx'^dx'^^Vi)} that arrises from (^, gives 

{vjVk ■ ■ - vi {dx„dx„Vi + dx'^dx'^Vi)} = {dx„dx„ + dx'^^dx'J (vj ■ ■ ■ v^) 

-2 {vk - ■ ■ vi [{dx„Vi) {dx„Vj) + {dx'^Vi) (dx'^Vj)] } , (37) 



where the right-most term in (^6|) has been subtracted from both sides of (|3g) to obtain (^7|). Note that 
{dx„Ui) {dx„Uj) = {dx„Vi) {dx„Vj) and {dx'^u'^) (dx'^u'j) = {dx'^Vi) (dx'^Vj), and that use of (||) gives 



{dx„dx„ + dx'^dx'J = 2 ( dr^dr^ + ^9x,.i9x, 



Then, (p7|) can be written as 



{vjVk ■ ■ - vi {dx^dx^Vi + dx'^dx'^Vi)] = 2 ^9r„9r„ + ^dx^dx^ {i 



Vi) 



-2 {vk--- VI [{dx^u,) (dx^uj) + (dx'yi) } . (38) 
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VI. APPENDIX B: DERIVATION OF (14) 



The purpose of this appendix is to derive (p^). Since allows Pi to be written as Pi — dxi {p — p'), the 
differentiation chain rule gives 

VjVk ■ ■ ■ viPi = VjVk ■ ■ ■ vidx, {p - p') = dx, [vjVk ■ ■ ■ vi {p - p')] - (p - p') {{dx,Vj) Vk ■ ■ ■ u/} , (39) 

where the notation {°} /i/ denotes the sum of all terms of a given type that produce symmetry under interchange of 
each pair of indexes with the index i excluded. Recall that the product VjVk • ■ ■ vi consists of iV — 1 factors. Sum 
the N equations of type ( ^ ) such that the sum is even under interchange of all pairs of indexes; then 

{v,Vk ■ ■ ■ viP,} = {dx^ [vjVk ■■■viip- p')]} -{N-l)ip- p') {{dx^Vj) Vk---vi}, (40) 

where use was made of the fact that the — 1 terms in the sum {{dxiVj)vk ■ ■ ■ '^i} /i/ each give the same result, 
namely, {{dxiVj) ■ ■ ■ vi}. From (^, dxiVj = dx^Uj — dx'.u[\ such that the definition of strain rate ( p^ gives 

{dx^v,+dx,v,) /2^s,,-s[,^. (41) 

Use of (^) gives {{dx^j)vk ■ ■ ■ vi} ^ {{{dx,Vj)vk •••«/} + {{dx,Vi) Vk ■ ■ ■ vi}) /2 = {(sy - Vk ■ ■ ■ w/}, substi- 
tution of which into (^) gives 

{vjVk ■■■viP,} = {dx^ [vjVk ■■■viip- p')]} -{N-l){p~ p') { (sy - s'y.) Vk---vi}. (42) 

Alternatively, consider that (|^) allows Pi to be written as Pi = 28^ {p + p')- Analogous to (|4^), it follows that 

{v,Vk ■ ■ ■ VlP^} = 2 {dr^ [VjVk ■■■Vi{p + p')]} -{N-l){p+ p') {i2dr,Vj) ffe • • • f/} . 

The 2dr^Vj can be replaced by, 2driVj — dx^Uj + dx'u[, such that 

{£2dr^Vj)vk ■ ■ - vi] = [{{2dr^Vj)vk •••«;}+ {{2dr^Vi) Vk ■ ■ ■ I'l}) /2 = {{sij + s'^^) Vk ■ ■ ■ w/}. Then, analogous to 
(El), it follows that 

{VjVk ■ ■ ■ VlPi} = 2{dri ['"j'^k ■ ■ - vi {p + p')]} - (iV - 1) {p+p') {{s^j + s'y) Vk - ■ - vi} . 



VII. APPENDIX C: THE COEFFICIENT IN (21) 

The purpose of this appendix is to obtain a formula for evaluation of the coefhcient {W[7v_2P]<5[2P] } in the 
special Cartesian coordinate system. In this coordinate system, ^ = Su such that W[jv] is the product of N Kronecker 
deltas of the form Sn. Consider setting A^i of the indexes equal to 1, an even number A'2 to 2, and an even number 
to 3 such that = A^i + iV2 + -^3. Then {W[jv-2P]'^[2P] } becomes a sum of zeros and ones. What is that sum? If 
all indexes are set to 1, i.e., Ni = N, then all terms in the sum {W[n-2P]S[2P]} ^-re unity such that {W[n-2P]^[2P]} 
equals its number of terms; from ( ^9| ) in Appendix E, that number is {x^2p) ^ -'^)''- "^^^ notation (^i^2p) ^ 
binomial coefficient. To interpret {2P — 1)!!, recall that g!! = g (g — 2) (g — 4) • • • ql, where ql is 2 or 1 for q even or 
odd, respectively, and ( — 1)!! = 1. Now consider setting two indexes to 2, thus iVi = — 2, and N2 = 2. Name the 
two indexes i = 2 and j — 2. The only term in {W[jv-2P]<5[2P] } that is nonzero is that which has i and j together 
in a single Kronecker delta, 5ij, within S^2P]- For P — there is no such 6ij, in which case {W[jv]'5[o]} = 0- For 
P > 1, there is one such 6ij, which is set to 1 and it multiplies the quantity {W[jv-2P]<^[2(p-i)] }; since this quantity 
is evaluated with all Is, it is equal to its number of terms, namely (''^ (^'.^pf "^'') (^ ~ 1) ~ l)-'- Now consider 
setting four indexes to 2; thus TVi = — 4, and A2 = 4. Name the four indexes i — 2, j — 2, k — 2, I ^ 2. The only 
terms in {W[7v-2P]'^[2P] } that are nonzero are those that have factors SijSki, SikSji, or SuSjk within S[2P]- For P < 1 
there is no such pair of Kronecker deltas such that {W[7v](5[o]} = and {W[^r_2]'5[2] } = 0. For P > 2, there are the 
above 3 — (A2 — 1)!! nonzero factors and each multiplies the quantity {W[jv-2P]'^[2(p-2)] }; since this latter quantity 
is subsequently evaluated with all Is, it is equal to its number of terms, namely (^^"^j^^p^"^'') (2 (-P — 2) — 1)!!. 

Continuation of this study reveals the pattern that when {W[^r_2P]<5[2P] } is evaluated with A^2 2s and A^i Is such 
that N ^ Ni+ N2 then 
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if 2P < N2 then {W[7v-2P]^[2P] } = 0, otherwise 
(iV - 2P) + 2 (F - ^ 



{W[jv-2Pi<5[2P]} = (A^2-l)!! 

= (iV2 - 1)!! 



(7V,_2(P-f))!(2(P-f))! 



2(P-^ I -1)!! 



If one ceases increasing the number of 2s and commences increasing the number of 3s in pairs 
such that N = iVi + iV2 + TVs, then (7V2 - (2 (P - ^) - l)!! is replaced by 

{N2 - 1)!! (iVa - l)!!('^"^'^Yi^*l^pP"^'') - ^ ^ 3^) ^ l)-- Of course, the binomial coefficient can be ex- 

pressed as follows: ('^"^'^^w-^pP"^'') = - ^^2 - N^Y-I [{N - 2P)! (2P - N2 ~ TV3)!]; also, (iV - iV2 - TV3)! = 
N\\. The double factorial can be eliminated by means of the following identities: 



(2Q-1)!!/(20)!-1/(2Q)!!-1/(2'3Q!) 



That is, (2P-iV2-TV3-l)!!/(2P-A^2-TV3)! = 1/ (2P - Ar2 - TV3)!! = 1/ 
{N2 - 1)!! = N2\/ [2^=/2 (7V2/2)!] . Finally, 

if 2P < + ^3 then {W[Ar_2P]'5[2P] } = 0, otherwise 

IV 2 2 



2(^-^-^)(P-f -f )! 



{W[^_2P]5[2P]} = (^2 - 1)!! (A^3 - 
= Ni\N2\Ns.\/ 



[N -2P)\2 
N2\. /Ns 



2 2 



(TV-2P).2n^l.(^l!fP-f -f ). 



(43) 
also, 

(44) 
(45) 



VIII. APPENDIX D: PROPERTIES OF ISOTROPIC SYMMETRIC TENSORS 



The purpose of this appendix is to determine (1) which components of a symmetric, isotropic tensor are 
zero; (2) how many components are zero and how many are nonzero; and (3) the relationships between the nonzero 
components. 

Consider which components -D[iV:A'i,JV,.iV3] are nonzero and which are zero. In (|2|), W[^_2P] vanishes if any 
of its indexes is 2 or 3, 5[2P] vanishes unless it contains an even number of indexes equal to 2 and a likewise even 
number of 3s (and of Is). Thus, a component i?[Af:Afi,Af2,W3] is nonzero only if both N2 and N3 are even. Because 
iV = A^i + iV2 + N3, -D[^.^j jv2,Af3] is nonzero only when iVi is odd if N is odd and only when A^i is even if N is 
even. The values of A^i that can give nonzero values of P'[Ar:Ari,Ar2,Af3] are N, N — 2, • • •, or 1; i.e., M + 1 cases 
of A^i. Given A^i, the values of [iV2, N3] that give nonzero values of D[N:Ni,N2,N3] are [N — iVi, 0], [N ~ Ni ~ 2, 2], 
• • • , [0, iV — TVi] ; i.e. , {N — Ni +2) /2 cases (note that N — Ni — N2 + N3 is necessarily even) . Counting the number of 
cases of [A^2, A^3] as Ni varies from to or 1, (i.e., substituting iVi = N, then iVi = TV - 2, • • • into (TV - iVi + 2) /2 
and adding the resultant numbers) shows that there are 1 + 2 + 3 + • • • + (A/ + 1) = {M + 1) {M + 2) /2 components 
D[N:Ni,N-2,N3] that are nonzero. Since there are 3^ components of D[jv], the remaining 3^ — (M + 1) (M + 2) /2 
components are zero. Since there are M + 1 linearly independent components (that are related to the Dpf p), there 
are {M + 1) {M + 2) /2 — {M + 1) = M {M + 1) /2 relationships among the nonzero components £)[A':iVi.iV2,Af3] • For 
instance, interchange of the values of N2 and A^3 produces components that are equal. 

Consider the M + 1 linearly independent equations that determine the Dn p in terms of M + 1 of the 
D[N:Ni,N2,N3]- Dpifi is related by (|2l| ) to only the component D[js[:n,o,o] because (^) shows that the coefficient 
of D]\[ Q, namely {W[jv]5[o]}, vanishes unless A^i = N. That is, if all indexes in (|2l| ) are 1, then D[js[:n,o,o] appears 
on the left-hand side and the D]\j^p for all P appear in the equation. This equation is essential for determining D^^q 
and is called "the equation for _Djv,o"; similar terminology "the equation for D^.p^^ is used below. With the equation 



for DNfi in hand, consider Dn,i- Dn,i is related by (21) to Z?[jv 



:Ar-2,2,0] 



or D 



[N:N -2,0,2]- 



When D 



[N:N-2,2,0] 



D[N:N-2,o,2] is on the left-hand side of ( |2l| ) the equation for D^,! results because the coefficients [n -2P]5[2P]} of 
Dn,p for P > 1 do not vanish, but the coefficient of Djq q does vanish. Now consider an equation for P'7v,2- Disi,2 



is related by ( |21| ) to i^[jV:Af-4,4,o] or I?[Ar:7v-4,o,4] or D\^is!-M-i,2,2\\ these components also involve Dj^^p for P > 3 but 
not for P < 1. This procedure repeats until the last equation is produced; only Dn^m appears in the last equation. 
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If N is even then D^^m is related to I)[Ar:o,Ar2,^3] with N2 and N3 equal to any positive even numbers such that 
N = N2 + N^. If N is odd then Dm,m is related to i^jjva.ATj^Ara] with N2 and equal to any positive even numbers 
such that N = 1 + N2 + N3. This procedure results in a set of M + 1 linearly independent equations that can be solved 
to obtain the Dm.p in terms of the D^fq-jq^ fq^ j^^y Note that M + 1 components must be chosen for use in the M + 1 
equations. For instance, from the above example, -D[Ar:Ar,o,o] must be used; either -D[Ar:JV-2,2,o] or ^[Af:Ar-2,o,2] must 
be chosen, and one of £'[Ar:jv-4,4,o] or -D[Ar:jv-4,o,4] or -D[iV:A'-4.2,2] must be chosen, etc. For simplicity, the chosen 
components can aU have A^3 = 0, i.e., the choice can be £'[Ar:Ar,o,o] , -C'[Af:Af-2,2,o]j -C)[jV:Ar-4,4,o] , • ' -D[JV:Ar-2M,2M,o]- 

The above procedure also reveals algebraic relationships between the nonzero -D[jV:Ari.7V2,Af3]- The equation 
for Djq^i can be expressed in terms of either -D[Ar:Ar_2,2.o] or -D[Ar:Ar_2,o.2] ! the left-hand side is the same in either 
case because the coefficients ( p3| ) are the same; hence -D[Ar:Ar-2,2,o] = -C'[7V:Ar-2 0,2]- The equation for Z?Ar,2 can be 
expressed in terms of -D[Ar:7v-4,2,2] or -D[Ar:Ar-4,4,o] or i?[jv:Ar_4,o,4] such that (|2^j gives -D[^.jv-4,4,o] = -C'[JV:W-4,o,4]; 
but what is the relationship of D^^.jq^i 2,2] to -D[Ar:7v-4,4.o] and D[N:N-i.0A]^ When £'[Ar:Ar_4,o,4] or I?[Ar:Ar-4.4,o] is 
on the left-hand side of @, the nonzero coefficients are, from (||), [N - 4)!4!0!/ [{N - 2P)!2^2!0! (P - 2 - 0)!] , but 
when D\^M:N-i,2,2\ is on the left-hand side, the nonzero coeflicients are {N - 4)!2!2!/ [{N - 2P)!2^l!l! (P - 1 - 1)!] . 
The ratio of these coefficients is 3. Since this ratio is independent of P, the entire right-hand side of is three 
times greater when -D[iV:Ar_4.o,4] is on the left-hand side as compared to when -D[Ar:Ar_4.2,2] is on the left-hand side. 
Therefore, the proportionality £'[iV:Ar-4.o,4] '■ -D[Af:Af-4,2.2] (and also P'[Ar:Ar_4_4.o] : -D[7V:Ar-4.2,2]) is 3 : 1. In general, 
for given TV and TVi, and hence given N2 + N3 = N — Ni, and another choice of iV2 and N3, call them N2 and iVg, 
such that N2 + = N2 + N3 = N — Ni, the proportionality obtained from ( p3|) is D^j^.j^^ j^^j^^] : D^j^.j^^ j^i jv'] = 
[A^zliVa!/ (iV2/2)! (7V3/2)!] : [iV^IA^g!/ (N^/^y. (^3/2)!]. Parameterized in terms of an integer L such that N = Ni + 2L, 

for given iVi, the proportionalities are P'[Ar:Arj,2L,o] : D[,m,,sj^^2L-2,2] ■ D[N:Ni,2L~iA] ■ ■■■ ■ D[N:Ni.o,2L] = [(2L)!0!/L!0!] : 
[{2L - 2)!2!/ (L - l)!l!] : [(2L - 4)!4!/ (L - 2)!2!] : • • • : [0! (2L)!/0!L!]. This constitutes {N - Ni) /2 relationships 
among the nonzero components D[N:Ni,N2.N3] for given iVi. Substituting the M + 1 cases of Ni (i.e., Ni = N, 
Ni = N — 1, • • •) into {N — Ni) /2 the number of relationships thus identified among the components of T^in] is 

+ 1 + 2H \-M = M (M -I- 1) /2. In the paragraph above, it was determined that the total number of relationships 

among the nonzero components of D[7v] is M (M 4- 1) /2. Consequently, all such relationships have now been found. 

IX. APPENDIX E: DERIVATIVES OF ISOTROPIC TENSORS 

The objective of this appendix is to develop succinct notation for isotropic tensors and their derivatives with 
specific attention to their first-order divergence and their Laplacian. Those derivatives appear in (p^). A derivation 
of those derivatives operating on an isotropic tensor that is symmetric under interchange of any pair of indexes is 
given. 

First, notation is developed: d[2p] is the product of P Kronecker deltas having 2P distinct indexes. For 
example, (5[g] = SijSkiSmn, where 6ij = 1 if i = j and 6ij = if i 7^ j. W[jv] is the product of N factors ^ each with 
an index distinct from the other indexes. For example, W[4] = ^-^^y- For convenience, define 

^[0] = 1, <5[_2] = 0, W[o] = 1, W[_i] = 0, W[_2] = 0. (46) 

The plural of 5 is 5s and that of W is Ws. It is understood that products of Ws (e.g., W[jv]Wj^]) and of 

i5s (e.g., S[2N]S[2K]) and of Ws with Ss (e.g., W[jv](5[2K]) have all distinct indexes. Then, the Ws factor, e.g., 

W[4] = W[i]W[3] = W[2]W[2] — W[i]W[i]Wp], etc., and the 6s likewise factor. The operation "contraction" means 

to set two indexes equal and sum over their range of values; the summation convention over repeated Roman indexes 
3 

is used, e.g., Su =^ Su — Su + 622 + 633 — 3. A contraction of two indexes of S[2P] produces either 3(5[2(p-i)] or 
1=1 

(5[2(p-i)] depending on whether the two indexes are on the same Kronecker delta or different ones, respectively. The 

2 

contraction of W[jv] on two indexes produces W[jv-i] because — = 1- Consider the contraction of ^ with 
W[jv](5[2P]- If the index i is in W[jv] then the contraction ^W[jv]<^[2P] is '^[n-i]6[2P] because = 1. If the index 

1 is in S[2P] then the contraction ^^[n]S[2P] is W[jv+i]'5[2(p_i)] because ^Sij = ^. The notation |[o|[^means the 
sum of N terms where each term contains a distinct index j, and the index j is interchanged with all implied indexes, 
but j is not interchanged with any explicit index. For example, 

$W[2,%$;^ = ^^5, + 'f'-lSa + 'fjS,,, (47) 

and Wril(5[2l -\\ = —Onm — — + —Snm 1 Ojm 1 <5„j , (48 



9 



where, i is an explicit index and is therefore not interchanged with j. 

French braces (i.e., {o}) means: add all such distinct terms required to make the tensor symmetric under 
interchange of any pair of indexes. For example, 

Note that terms that are necessarily equal do not appear; i.e., since ^^5^1 appears, neither ^^Ski nor ^^Sik 
appear. Because of the commutative law of addition, {o} commutes with addition; e.g., {Wj^v]} + {W[q]5[2p]} = 
|W[^] + W[Q](5[2p] }. Because of the distributive law of multiplication, multiplication by a scalar function commutes 
with the {o} notation; i.e., A{r) {W[7v](5[2p] } = {A {r)'W[j^]S[2p]} ■ 

The number of terms in various sums, {o}, is required repeatedly: {^[2P]} has (2P— 1)!! = 
(2P - 1) (2P - 3) (2P - 5) • • • (1) terms. Since W[jv](5[2P] has 2P + N indexes, {Wj^rjJpp]} has (^^^^) (2P - 1)!! 
terms, where the binomial coefficient (^^^'^) is the number of ways of selecting the N indexes in W[7v] from the 
total 2P + N indexes. If i is an index in {W[jv](5[2P] } , then {W[^](5[2P]} has (^^^^^) (2P — 1)!! terms in which i 
appears in a Kronecker delta because there are N indexes to select for Wj^v] from the remaining 2P + N — 1 indexes. 
Similarly, {W[n]S[2p]} has (2P- 1)!! terms in which i appears in a factor (r^ / r) because there are — 1 

indexes remaining to select for Wj^vj from the remaining 2P + iV — 1 indexes. Note that {Wj^vj} has only 1 term. 
Hence, {W[^]} = W[^]. In summary, 

/op + N\ 

{W[w]^[2P]} has f ^ j (2P-1)!! terms; (49) 

{W[jv]^[2P]} has (^^^^ ' ^) (2P - 1)!! terms with i in 6[2p]; (50) 

{W[jv]^[2P]} has (^^^ ^ ~ (2P - 1)!! terms with i in W[jv]. (51) 

The sum of the number of terms in (|5^|l]), namely f ^+/~^) (2P - 1)!! + (^-^+^,7^) (2P - 1)!! = 
(2P + TV _ 1)1 + j^j^ (2P - 1)!! = (^^^^) (2P - 1)!!, agrees with the total number of terms in 

Now, rules for differentiation of symmetric, isotropic tensors are developed. Note the identity 

dr. (rj/r) = [% - (r,r,/r2)] /r = (^p, - Wp]) /r, (52) 

from which it follows that 

dr, (ri/r) = (3 - 1) /r = 2/r, and {vj/r) = [r^ - (r^r^/r^)] /r = 0. (53) 

The latter formula greatly simplifies the divergence of W[jv] because operating on W^jy] vanishes when it operates 
on any factor other than the factor ^ within W[^]. 

The divergence and gradient of W^yv] are needed. If i is an index in Wj^v] then the divergence of Wj^v] is denoted 
by Vr • W[jv] = 9riW[jv]. Application of ( |5^ ) gives dn^yi] — {n/r) — 2/r, and d^'W^^] — W[7v_i]9ri (fi/r) = 
^W[7v-i]. The gradient of W^^r] is denoted by dn^iN] where i is not an index in W[jv]. From the differ- 
entiation chain rule, dr^'W^^-^ is the sum of N terms, each of which has the form W[jv_i]9ri {fj/f)- Therefore, 

by use of (||), a^.Wf^r] = $W[Ar_i] [5ij - {r,rj/r'^)] /r-^^^ = ]:'P^[N-i]5ij1l^ - ^^[n+i]i where use was made 
of $W[^._l] [-{nr,lr^)]lrf^ = [- (r./r) /r] $W[^._i] (r,/r)$f = - (r,/r) /r (iVW[^j) = -fW[^.+i], because 
^W[jv-i] (''j/^)$^ is the sum of N identical terms each equal to W[7v]. In summary, 

2 

If i is in W[7v] then Vr • Wjyv] = 9,.,W[jv] -W[Ar_i]. (54) 
If i is not in W[Ar] then 9r.W[jv] = - ^W^^^i^Sij^^ - — Wf^v+i], (55) 

Consider the divergence of ^[n]&[2P]- If the index i is in W^jv] then, from (]5^), 9,.. iyV[N]^[2P]) = 
(5[2P] W[jv|i] = |<^[2P]W[7v_i]. If the index i is in 5[2P], then (given that index k is not in W[jv]) {y^[N\5[2P]) = 
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SikdriW[N] = <5[2(P-l)]5rj^W[7v] = <5[2(P-l)]:F$W[Ar_i]4j:^ 



AT 



follows from (|55|). In summary, 
If i is in W[7v] then 9, 
If i is in (5[2p] then d, 



S[2{P~i)]^'^[N+i], where the last expression 



{^[N]Sl2P]] 



)] - P^[N-l]Skjt'j' - '^[2(P-1)] — W[Ar+i] 



[2(P-1 



(56) 
(57) 



The above results allow evaluation of the divergence {W[7v]^[2P] } = Vr • { W[7v](5[2P] }• It follows from use 
of (^6|) and the distributive law of multiplication and the fact that the number of terms in {W[jv]^[2P]} in which i 
appears in the factor ^ is the same as the number of terms in {W[n-i]S[2P]} [see ( p9[ - |5T| )], that for those terms in 
which i is in W[7v], the divergence of {W[7v](5[2p] } yields ^ {W[jv-i]^[2P] }• Similar use of ( p7| ) gives that for those 

terms in which i is in S[2P], the divergence of {W[jv](5[2P] } yields ^ {W[n-i]S[2P]} - ^^^'''^^ {W[jv+i]<5[2(p-i)] }• 
Thus, 



Vr • {W[Ar](5[2P]} = - { W[^_i](5[2P] } + 



— {W[N-l]Sl2P]} {W[^^+l](5[2(p_l)]} 



r 



1){W 



[Af-l]0[2P] 



} 



N{N +1) 



{W[7v+i](5[2(p~i)]} ■ 



(58) 



Because of the definitions in ([46[), (58) remains valid if TV is or 1 or if P is or 1. 

Derivation of the formula for the divergence of an isotropic tensor requires evaluation of the contraction 

(2P — 1)!! occurrences of the index i within 6[2P] and 



f2P+N-l\ 
N 



{W[N]6[2P]} 7". From (|50|), in {W[jv](5[2P] } there are 
each gives the contraction Sij ^ — which decreases P by unity and increases N by unity thereby producing several 

(2(P— 1) — 1)!! terms in a {W[jv+i]5[2(p-i)] }; thus the 



{W[^r+i](5[2(p_i)]}. From (|49|), there are ( 



.2(P-l) + (Ar+l)\ 



number of {W[jv+i](5[2(p-i)] } so produced is ( 



(JV+l) 
2P+A'-l\ 
N ) 



(N + l). 

From (|5l|), the contraction {WjjvjJpp] } ^ contains C^^-^^) (2^-1)!! ter ms in which i appears within Wj^v] and 



(2P-l)!!j/[f(^-;^+(f+i)) (2(P-1)-1)!! 

2P+7V-1\ 
N~l ) 

each results in the contraction = 1, which decreases N by unity. The number of {^[n-i]5[2P]} so produced is 
"(2^+^-1) (2P- 1)!!] / (2P- 1)!!] = 1 because {MV{n-i]S[2P]] has ('^^^i]'^) (2P- 1)!! terms, which 

is also the number of terms given in (|5l|) 



Thus, 



contraction on i: {W[n^5[2p]] — = [N + I) {W[Ar+i]^[2(p_i)] } + { W[Ar_i](5[2P] } 



(59) 



The general isotropic formula for a tensor ^[n\ (r) of order N that is symmetric under interchange of any pair 
of indexes is 



A[jv] (r) = Ao (r) {V^[n\] + (r) {^[N-2Ai\] + M (r) {W[Ar_4]'^[4] } + 



T, 



last 



M 

E 

P=0 



Ap (r) {W 



[W-2P]0[2P] } I 



(60) 



where the ^o(^), ('')i etc., are scalar functions of r, and Tiast is the last term. Note that for brevity in this 
appendix, the subscript N has been omitted from A^v.o {r), [r), etc. If N is even, then Tiast = ^7v/2 ('') {^[w]} 
and M = N/2. If N is odd, then Tiast = ^(Ar-i)/2 {r) { W[i]5[w_i] } and M = (iV - 1) /2. 

All of the foregoing has set the stage for efficient derivation of a formula for the divergence Vr • ^[n\ (r). Also 
needed is the fact that the gradient of a scalar function of r = y/r^rl is driA{r) = (dnr) drA{r) — ^drA{r) — 
Wt^drA ( r). Consider the divergence of a term in (|60|). By use of the differentiation chain rule, and substitution of 
( p9| ) and (jSq), 



Vr • [Ap (r) {W 



[Ar-2P]'5[2P]}] — 

{W[Ar_2P]^[2P] } -^drAp (r) + Ap (r) Vr • {W[7V-2P]'5[2P] } 
= [{N - 2P + 1) {W[Ar_2P+l]^[2(P-l)] } + {W[Ar_2P-l]^[2P] }] OrAp (r) 



+Ap (r) 



(P+l){W[Ar_2P-l]'5[2P]} 
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'Ap (r) 



(iV-2P)(iV-2P+l) 



2P+1]0[2(P 



-1)]} 



= Bn.p (r) {WjAT^ap+il^pcp-i)]} + Cp (r) { W[Ar_2P-i]^[2P] } , (61) 
where Bat^p (r) and Cp (r) are defined by the foUowmg operators, Ob {N,P) and Oc (P), operatmg on Ap (r): 

N~2P] 



Bn,p (r) = Ob {N, P) Ap (r) , where Ob [N, P) = {N ~ 2P + 1) 



dr 



Cp (r) = Oc (P) Ap (r) where Oc {P) 



Thereby, the divergence of (|6C|) is 

M 



dr + -{P+l) 

r 



M 



Vr • AfAT] (r) =^ Bn.p (r) { W[7v-2P+i]^[2(p-i)] } + ^ ^p (r) {W[Ar_2P-i]'^[2P] } 



(62) 
(63) 

(64) 



p=o 



P=0 



where (|l]) and (g|) were used. 

Now, (64) can be checked by comparison with the divergence performed on the explicit-index formulas for 
symmetric, isotropic tensors of rank 1 to 4. The lowest-order tensor for which the divergence is defined is a vector 
(i.e., N = 1), in which case (H) gives 



dr^A, (r) -a.Ao {r) + -Ao (r) , 
r 

which is easily verified by evaluating the divergence of a isotropic vector, namely [Aq (r) ^] . Expressed with 
explicit indexes as well as in the implicit-index form of (^0|) , isotropic tensors of rank 2 to 4 that are symmetric under 
interchange of any pair of indexes are: 



A,j (r) = Ao (r) -f + Ai (r) % = A^ (r) Wp, + Ai (r) <5[2] . 



(65) 



(r) = Ao (r) + (r) (^<5,fc + ^S,k + y<5.,) = (r) Wp, + A, (r) {W[i]<5[2]} 



(66) 



A.,ki (r) = Ao (r) ^-^^^ + A, (r) ( + '^5,i + '-^5^ + '-^S^u + '-^5,k + '-^5,, 



+A2 (r) (%4; -f- SikSji + SjkSii) 

Ao (r) W[4] + A, (r) { W[2]5[2] } + A2 (r) {6[^ } . 



7" 



TjTl 



(67) 



One can see the brevity of the implicit-index formula as the rank of the tensor increases. The first-order divergences 
obtained by differentiating the above explicit-index formulas as well as from (p3) are: 



Vr • A[2] (r) = 



dr + -] An (r) + drAi (r) 



[i?2,i(r) + Co(r)]{W[i]<5[o]} 



(68) 
(69) 



Vr • A[3] (r) = 



dr + ^ j [r) + [2dr - ^ Ml (r) 



dr + -] (r) 
r ' 



[B3,i (r) + Co (r)] {W[2](5[o] } + (r) { Wpj^^p, } 



(70) 
(71) 



Vr • A[4] (r) 



dr + ^] Ao (r) + (^39. - ^ Ml (r) 



rjrkTi 



dr + -]Ai (r) + 9,^2 (r) 



= [S4.1 (r) + Co (r)] {W[3]<5[o] } + (i?4,2 (r) + Ci (r)) {W[i]<5[2] } 



(72) 
(73) 
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In the implicit-index formulas in (|69|; |7l|, [73|), terms from ( |64| ) that are zero because of (|4^) have been omitted. 
Equation (H) has been checked by using the implicit-index formulas in ( |69| , 71 , |73| ) to obtain the explicit-index 
formulas in (|68|, |70| , fz^ ), respectively. 

The Laplacian of a symmetric, isotropic tensor is also needed for the term 2j/VrD[^] in (|l^). Application of 
@ to V2W[jv_2P] and use of (|5|-||) gives 



[N-2P] = {W[^r_2P-l]a,.„9,.„^} +2{W[^r_2P-2] (^-.y) (^--.7)} 



(74) 



Now, {W[7v-2P-i]9r„9r.„^} is ( ) terms, each one is of the form W[^r_2P_l]5r„9,.„ ^ = W[^r_2P-i] (t^t! 
!-W[7v_2P]- Also, {W[^_2P-2] (<9r„^) '7") } IS (^^^"^) tcrms, each one is of the form 



W[jv_2P-2] (5.„^) (9,.„^) = W[^,_2P-2];^(4,- 
{'W[Ar_2P-2]<5[2]} has (^^^^) tcrms; thus, @ gives 



;^W[Ar_2P-2]% - ;:^W[^r_2P]. From (|4|), 



V2W[Ar_2P] = -J {W[Ar_2P-2]'5[2]} " -J 



N - 2P\ fN - 2P 



W 



[W-2P] 



(75) 



The binomial coefficients prevent a nonzero term in ( [Tq ) when W[7v_2P-i] or W[7v_2P-2] vanish in (Q) as required 
by definition mm provided that we define 



N ~2P 



if iV- 2P < 1, and 



N - 2P 



if TV - 2P < 2. 



(76) 



Of course, (^6|) is consistent with 1/K\ = for K < (Abramowitz and Stegun, 1964, equation 6.1.7). Given ([7^), 
we can define, for brevity 



N-2P 



N ~2P 



N ~2P 
1 



Now (111) and (0) give 



V^W[^_2P] — — {W[7V-2P-2]<^[2]} - ^ 



2 W[^_2P]- 



(77) 



(78) 



Now, use of (|^) gives 



2 {^lN-2P-2]d[2]} 



[Ar-2P] 



5[2P] 



- -2 {{'^[N-2P-2]S[2]} S[2P]} -2—{'^[N-2P]S[2P]} 



r 

R{N,P) 



(79) 
(80) 

Noting the appearance of W[^_2P-2] in ([^ and recalling that W[jv-2P-2] = if — 2P — 2 < 0, the coefficient 
R {N, P) is defined by 



2 { W[Ar_2(P+l)]'5[2(P+l)] } - "^^^2^^ { W[Ar„2P] ^[2P] } 



R{N,P) = if A^ - 2F - 2 < 0, 



otherwise, R [N, P) =2 



N-2P 
2 



( " 

\2P + 2 



(2P+ 1)! 



(81) 
(82) 



The coefhcient (^^^0/ {2P+2) (2P+ 1)!! is the number of terms in {W[jv -2P-2]'5[2]} divided by the number in 
{W[jv-2(p+i)]'5[2(p+i)] } as obtained from (4^. Because of (p3|), (81-^) can be simplified to 



R {N, P) = [2^+1 (A^ - 2P)\ {P + 1)!/A^!] 6 (A^ - 2P - 2) , 
where Q {x) — 1 for a; > and Q [x) — Q for a; < 0. 



(83) 
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The Laplacian of the product of two functions / and g is (^3[). When appHed to (|60|), the case f = Ap (r) 
and g = {W[7v-2P]^[2P] } are needed. Recall that driA{r) = ^drA{r). The last term in (|33D vanishes as follows: 
{drj){dr,g) = [IdrAp {r)\ridri {^[N-2P'fi]2P]) = {[TidTyf[N~2P]\5[2P]} 0; this vanishes because 

(H) shows that r,d.ryf[N-2P\ = 0. Then (|o|) used in (||) combined with V^A (r) ^ [df. + f 9^) A (r) give 



V2 [^p (r) {W[N-2P]S[2P]}] 



dl + -a, 

+Ap (r) 



2P 



R{N,P) 



Ap [r) 



{W[7V-2P]<5[2P]} 



{W[Ar_2(P+l)]^[2(P+l)]} 



(84) 



M 



The Laplacian operation on ( |60| ) is simply the sum, ^ , of terms (p4[). 

p=o 



X. APPENDIX F: MATRIX ALGORITHMS 



For computations, it is useful to write (21) as a matrix equation. Let the column index be J = P + 1, and 
the row index be / = {N2/2) + 1, such that both J and / range from 1 to M + 1 in (|2^). Use 7V3 = in (|2[]2|) to 
define the following matrix elements 



Mm {I, J) — 0, for J < I, i.e., (I, J) = below the main diagonal; 

Mn {I, J) ={N-2I + 2)1 {21 - 2)1/ [{N - 2 (J - 1))12''-^ {I - 1)1 (J - J)!] , for J > /. 



(85) 
(86) 



The chosen M + 1 linearly independent components of D[jv] are arranged in a column vector having -D[Ar:Ar_2/+2.2/-2,o] 
as in its I-th row, and the M + 1 scalar functions -Dat.p are likewise arranged in a column vector having ZJjv./-! in 
its I-th. row. Then ([2l|) is written as the matrix equation 



/ D^N-.NfiM \ /Mjv(l,l) Mjv(l,2) 



^[Ar:Ar-2,2,0] 
\ -D[Ar:Ar_2Af,2M,0] / 



Mn (2, 2) 



Mjv(l,M+l) \ 
Mn (2,M + r 



V 





/ Dn.o \ 




Dn.i 


/ 


\ Dn,m J 



(87) 



••• Mjv(M + l,AfH 

Denote a matrix having matrix elements A (/, J) by 1^ i^i J)\ Then (|7|) and its solution are (respectively) 

|-^[Af:Af-2J+2,2J-2,0]| = lAfjV {I , Tj\\L> N ,1 ~l[ a,nd \^NJ^ = \Mn (/, J)r ^ P[jV:jV-2J+2, 2.7-2, 0]| . 

where \Mn (I, JT ' is the inverse 

of |Mjv (/,/]! 

The determinant of 

IMjv (/, J)| 

is the product of its diagonal elements; 

from ( |8^ ) that product is nonzero, hence exists. This inverse matrix is to be calculated numerically. In 

effect, evaluation of the components D^N:N-2j=t2,2J-2,o] by means of experimental data or DNS data and use of the 
solution in (^) produces the Dn,p for use in (|2ll) to completely specify D[jv]- 

A matrix algorithm is useful for determining the isotropic formula for the first-order divergence Vr • D[7v+i]. 
By replacing by + 1 and the symbol ^ by D in the divergence formula (|6^ ) , we have 



M' 



Vr D 



[JV+l] 



^ {W[^_2(p_i)](5[2(p-i)]} Oi3 {N + 1,P)Dn+i,p 



P=0 
M' 



Ob (iV + l,P) 



^ {W[jv_2P]'^[2P]} Oc (P)P'Ar+l,P, 

(A + 1) - 2P 



p=o 

((A+1)-2P+1) [dr 



Oc (P) 

M' = N/2 if N is even, and M' = 1 + (A - 1) /2 if A is odd. 



(89) 

(90) 

(91) 
(92) 
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The differential operators, i.e., dr = 9/9r, in (|9C|-|9l[) are obtained from (|62|-|63|) , and ( |92| ) is obtained by replacing 
iV by + 1 in ( po|) and simplifying and rearranging the terms. Comparison of (P2| ) with (20) shows that if N is 
even then M' — M; thus the matrix representation of { W[7v_2P]^[2P] } within (|8g|) is the same as in (p7|), which 
representation was abbreviated by \Mn(I,J)\ above. On the other hand, if N is odd, then M' = M + 1, and the 
last column of the matrix representation of {^[n-2P]^12P] } within ( p9| ) corresponds to P = M' = 1 + {N — 1) /2, in 
which case {'W[n-2P]^[2P]} contains Wj_]j== such that the last column of the matrix is zero. Thus, the matrix 
representation of { W[7v_2P]'5[2P] } within ( |89| ) is 



( ° 

\Mn {I, J)\ \ I if is odd; 

V 



if TV is 



(93) 



(94) 



In addition to the coefficient {W[7v-2P]<^[2P] } , ( p9| ) contains the coefficient {W[jv_2(p-i)]<J[2(p-i)] }■ From 
the matrix representation of {W[jv-2P]'5[2P] }? namely (|8q^p6|), the matrix representation of { W[7v-2(p-i)]'^[2(p-i)] } 
is (recall that J = P + 1) 



M'j^ (/, J) = 0, for J — 1 < /, i.e., (/, J) = on and below the main diagonal; 

M'j^{I,J) = (iV-2/ + 2)!(2/-2)!/ [(iV - 2 (J - 2))!2'^-2 (/ - 1)! (J - 1 - /)!] , for J - 1 > /. 



(95) 
(96) 



The matrix having these elements is denoted by Because of (^), if N is odd, then the matrix 

contains the matrix \Mn (I, J)\ shifted to the right by one column and a first column of zeros is included; that is. 



\M'^{I,J)\ = 



/O 



|MAr(/,J-l)| 



if N is odd. 



(97) 



Because of (^), the same is true if A'' is even except that the right-most column of \Mn (/, J)| is discarded. Thus, 



Ennil= I : | M^,(/,J-1) | 




if N is even (discard the right-most column). 



Define operator matrices that are of dimension M' -I- 1 by M' + 1, that have zeros off of the diagonal, and that 
have the operators ( p0[]9l| ) on the diagonals. Thus, recall that J = P + 1, and that dr = d/dr, and define matrix 
elements 



B (/, J) = Si J iN-2J + 4)[dr- 



N-2J + 3 



and C (/, J) = 6ij [dr + — 



(99) 



The matrices corresponding to Ob [N + 1, P) and Oc (P) in (|90|-^l|) are denoted by 1-5 (/, T^ , and|^lZiZ3, respectively. 

Let the components of Vr • D[jv+i] be denoted by (Vr • D[jv+i]) j^.j^^ ^^j, 
Vr • D[^_|_i] is a tensor of order N. In matrix notation, ( ^9| ) gives 



which denotes the fact that 



/ i^--'^lN+l])[N:N.,0,Q] \ 
(Vr • D[Ar+i])j^^^_2,2,0] 



V (Vr -DlAT+i]) 



(/, j)i mT;ji\ + mN^,j)\ icTTJ]]" 



[N:N-2M,2M 



a) 



Dn+is 



\ Dn+i,m' J 



(100) 



When apphed to D[Ar+i], the solution of (0) is Pn+ij~i\ = \Mn+i {I, J)] ^ |-^[jv+i:W+i-2.7+2,2,7-2,ol ^ substitution 
of which into (|100D gives 
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\ (Vr • D[Ar+i]' 



[Af:Af-2M,2M,0] 



/ -D[7v+l:W+i,o,0] 
^[7V+l:W+l-2,2,0] 

\ -D[7V+l:Ar+l-2M',2Af',0] / 



where, EZ23^ [1^;^ (/, j)ipT7;j^ + |Af^ (/, j)iPTr;7^] IMat+i (/, j)v \ 



(101) 



(102) 



We see that is the operator matrix that operates on the column matrix representation of D[jv+i] to produce 

Vr • D[7v+i]; this is true for any completely symmetric isotropic tensor, not just true for Dj^v+i]- 

It is helpful to illustrate this algorithm for N ^ 2 and = 3. Two examples are needed because the algorithm 
differs for even N as compared to odd N. For N — 2, (102) is 



Af2(l,l) 

Af2(l,l) Ma (1,2 
Af2(2,2 



5(1,1) 
5(2,2) 
C(l,l) 
C(2,2) 



M3(l,l) Af3(l,2) 
Ms (2, 2) 



dr + ^ 

' r r 
dr + t 



Computer evaluation of (|103|) produced (104). Consequently, (101) is 



(^r •D[3])p^2,0,0] 
(^r-D[3])p^g2,0] 



a, + 2 -4 

dr + ^. 



[3:3,0,0] 
[3:1,2,0] 



{dr 



7) Dill — 7-D122 



122 



where explicit-index notation is given at far right by use of £'[3:3,0.0] = ^111 a-nd i?[3:i.2,o] = ^122- 
For iV = 3, \Y {I, J)\ from (p^ is 



(103) 
(104) 

(105) 



Af3(l,l) Af3(l,2 
A/3 (2, 2) 



Ms (1,1) Ms {1,2) 
Af3(2,2) 



dr + l -f 

dr + t-^. 




5(1,1) 
5(2,2) 
5(3,3) 

C(l,l) 
C(2,2) 
C(3,3) 



Af4(l,l) A//4(l,2) A//4(l,3) 
Af4(2,2) Af4(2,3) 
M4 (3, 3) 



As with (104), the matrix was evaluated using a computer program. Consequently, ( |lOl|) is 



(^r-D[4])j3^3jj,jj 



Vr-D 



[4] j [3:1,2,0] 



I J. J. 

dr + 





4 _ 
r 3r 



^1 

Di 




_( {dr 
{dr 



'^^ Dun — 7^1122 



i^ii 



^-D2222 



(106) 

(107) 



where explicit-index notation is used in ( |107| ). 

A matrix algorithm is also needed for the Laplacian of a symmetric tensor. Performing the Laplacian of ( |60| ) 
and use of (^ gives 



[N] 



M 



{W[N-2P]S[2P]} {df.+^dr~ 



Sn-2 



P=0 \ + {W[jv_2(p+i)](5[2(P+l)]} ^'"'^'i^^ Dn,P 



D 



N,P 



(108) 



It is necessary to recall the definitions (Q) and 
is the same as in (B5- 



obtained from (pa-36) by replacing J by J + 1, i.e. 



The matrix representation of {W[n-2P]S[2P] } within ( |lO^ ) 
namely \Mn {I, J)\ . The matrix representation of { W[jv-2(p+i)]^[2(p+i)] } within ( |108| ) is 
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M* (/, J) = 0, for J + 1< /, 

M* {I, J) = {N -2I + 2)\ (2/- 2)!/ [(A^ - 2 J)!2'^ (/ - 1)! ( J + 1 - /)!] , for J + 1 > /. 



(109) 
(110) 



This is just the square matrix that appears in (^^ except that the left-most column in ( p7| ) is discarded and the 
matrix is then shifted leftward by one column and the right-most column is zeros. Those zeros appear because in 
the right-most column J = P+l^M + l such that {I, M + 1) contains the factor 1/ (TV - 2 {M + 1))! which is 
1/ (—2)! = if is even and is 1/ (—1)! = if A^ is odd(see Abramowitz and Stegun, 1964, equation 6.1.7). Thus, 

/ 0^ 

(discard the left-most column). 



M^(/,J) 



Mn{I,J+1) 



\ 







and 



Define two operator matrices that are zero off the main diagonal and contain ydf + ^dr ~ ^"^2^ 

on the main diagonal; i.e., their matrix elements are E{I,J) — Sij (j)^ + — ^ and F{I,J) 

6ijR {N, J — 1) /r^- Analogous to the derivation of (|10lD, the matrix representation of (|108D is 



R{N,P} 



V ^l'^m)[N:N~2M,2M,0\ J 

where, 



( D 



[N:N,0,0] 



\ 



[N:N-2,2,0] 
\ D[N;N^2M,2M,0] / 



\mn {I, j)\ wit:j)\ + Wifji, J)\ nryji 



\mn (/, j)r 



For both N — 2 and A = 3, the matrix representation of IS 



For A^ = 2 (111) is 



A/^(l,l) Mat (1,2) 
Mat (2, 2) 
Mjv(l,2) 
Afjv(2,2) 



E{1,1) 
£^(2,2) 
£^(1,1) 
F(2,2) 



Mjv(l,l) M^(l,2) 
Mat (2, 2) 



(^?I5[2]) [2:2,0,0] 
(V?D[2])[2^„2_pj 



2 a 4 



7^ 

2 Q 2 



D22 



mm 



D[2:2fl,Q\ 
^[2:0,2,0] 



■dr - 4) + ^D22 



(V?D[2])n 

(Vr'D[2])22 



22 



where the matrix was evaluated using a computer program. For A^ = 3 the matrix algorithm is 



(V?D[3]),22 



Dm 
-D122 



'■2 
7^ 



8 



12 
-pi 

dl 



Dm 
D122 



72- -Dill + ( 



-Dill + ^Di22 



^dr)D. 



(Ill) 



(112) 



(113) 
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